#! /usr/bin/env python
# Example scipt to show integration of a 1D spectrum

import nmrglue as ng
import numpy as np
import matplotlib.pyplot as plt

# read in the data from a NMRPipe file
dic,data = ng.pipe.read("1d_data.ft")
length = data.shape[0]

# read in the integration limits
peak_list = np.recfromtxt("limits.in")

# determind the ppm scale
uc = ng.pipe.make_uc(dic,data)
ppm_scale = np.linspace(uc.ppm(0),uc.ppm(data.size-1),data.size)

# plot the spectrum
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ppm_scale,data,'k-')

# prepare the output file
f = open("area.out",'w')
f.write("#Name\tStart\tStop\tArea\n")

# loop over the integration limits
for name,start,end in peak_list:
    min = uc(start,"ppm")
    max = uc(end,"ppm")
    if min>max:
        min,max = max,min

    # extract the peak
    peak = data[min:max+1]
    peak_scale = ppm_scale[min:max+1]

    # plot the integration lines, limits and name of peaks
    ax.plot(peak_scale,peak.cumsum()/100.+peak.max(),'g-')
    ax.plot(peak_scale,[0]*len(peak_scale),'r-')
    ax.text(peak_scale[0],0.5*peak.sum()/100.+peak.max(),name,fontsize=8)

    # write out the integration info
    tup = ( name,peak_scale[0],peak_scale[-1],peak.sum() )
    f.write("%s\t%.3f\t%.3f\t%E\n"%tup)

# close the output file and save the plot
f.close()
ax.set_xlim(ppm_scale[0],ppm_scale[-1])
fig.savefig("plot.png") # change this to plot.pdf, plot.ps, etc
